A meta-analysis fo the multi-trophic effects of grazing

Grazing effects are multi-trophic through mediating plant composition: a meta-analysis

Alessandro Filazzola, Batbaatar Amgaa, Charlotte Brown, Issac Heida, Jessica Grenke, Margarete Dettlaff, Tan Bao, & JC Cahill

library(tidyverse)
library(PRISMAstatement)

Purpose

Conduct a meta-analysis of the literature testing the indirect effects of grazing on animal taxa’s through the direct effects on the plant community.

Objectives

Timeline

date	task
Nov 9	Establish search terms to be used in the meta-analysis
Nov 12	Compile list of journal artcles and sub-divide for each researcher
Nov 14	Begin reviewing papers and extracting data
Jan 28	Complete data extraction from papers
Feb 11	Complete preliminary analysis and set structure for MS
Feb 25	Settle on analyses to be used and begin writing manuscript
March 11	Complete first draft of MS and pass to co-authors
March 25	Comments passed back on draft
April 2	Complete revisions and submit to journal

Literature Review - 1. Search

Search Terms

A systematic literature search was conducted using Web of Science for all emperical research articles. The review will include all studies globally. The intended purpose of this search is to capture all articles that have documented grazing either along a gradient (e.g. different frequencies or intensity) or presence/absence (e.g. excluded, ungrazed, or retired ranch lands). We condcuted two separate searches to capture studies that tested gradients and studies that compared grazing to ungrazed treatments. Duplicate articles between the searches were removed. We also intentionally excluded terms that resulted in articles not relevant to the purpose of this study including: review, synthesis, policy, social, carbon, and fish. The search terms that used were:

Search A graz* OR livestock AND exclosure* OR exclusion OR exclude* OR ungrazed OR retire* OR fallow*

Search B grazing intensity OR grazing gradient OR stocking rate OR rotation* grazing

Literature Review - 2. Sort

This steps includes a. checking for duplicating, b. reviewing each instance for relevancy, c. consistently identifying and documenting exclusion criteria. Outcomes include a list of publications to be used for synthesis, a library of pdfs, and a PRISMA report to ensure the worflow is transparent and reproducible. Papers were excluded with the following characteristics:

• Not emperical study (e.g. review, book chapter)
• Irrelevant categories (e.g. political science, law, sports tourism, art)

Reasons for exclusion

evidence <- read.csv("data//synthesisdata//evidence.csv")
### Identify studies that were excluded
excludes <- evidence %>% group_by(reason.simplified) %>% count(exclude) %>% filter(reason.simplified!="")
ggplot(excludes, aes(x=reason.simplified, y=n)) + geom_bar(stat="identity") + coord_flip()

## frequency of study
year.rate <- evidence %>% group_by(Publication.Year) %>% summarize(n=length(Publication.Year))

ggplot(tail(year.rate,30)) + geom_bar(aes(x=Publication.Year, y=n), stat="identity") + ylab("number of published studies") +xlab("year published") +theme(text = element_text(size=16))

Prisma report

## total number of papers found
nrow(evidence)

## [1] 2989

## number of papers found outside of WoS
other <- read.csv("data/synthesisdata//other.sources.csv")
nrow(other)

## [1] 0

## number of articles excluded
excludes <- evidence %>% filter(exclude=="yes")
nrow(excludes)

## [1] 2679

## relevant papers
review <- evidence %>% filter(exclude!="yes")
nrow(review)

## [1] 310

## papers for meta
datasets <- read.csv("data//binary.simplified.csv")
meta <- length(unique(datasets$uniqueID))
meta

## [1] 216

prisma(found = 2989,
       found_other = 1,
       no_dupes = 2989,
       screened = 2989,
       screen_exclusions = 2675,
       full_text = 315,
       full_text_exclusions = 0,
       qualitative = 315, 
       quantitative = 216,
       width = 800, height = 800)

## Warning in prisma(found = 2989, found_other = 1, no_dupes = 2989, screened
## = 2989, : After screening exclusions, a different number of remaining full-
## text articles is stated.

Patterns of GI Studies Globally

require(ggmap)
###  Start with base map of world
mp <- NULL
mapWorld <- borders("world", colour="gray50", fill="gray50") # create a layer of borders
mp <- ggplot() +   mapWorld

## colorblind-friendly palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#000000")

meta <- read.csv("data//synthesisdata//evidence.csv")
meta <- meta %>% filter(lat != "") ## remove blanks for lat and lon
meta <- meta[!grepl(";", meta$lat),] ## remove multiple entries
meta$lat <- as.numeric(levels(meta$lat))[meta$lat] ## convert from factor to number
meta$lon <- as.numeric(levels(meta$lon))[meta$lon] ## convert from factor to number


## plot points on top
mp <- mp+ geom_point(data=meta , aes(x=lon, y=lat), size=2) 
mp

3. Synthesis - Abundance

meta <- read.csv("data//binary.simplified.csv")
  
  
## Load packages and functions
library(reshape2)

## Warning: package 'reshape2' was built under R version 3.5.2

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

library(metafor)

## Warning: package 'metafor' was built under R version 3.5.2

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:tidyr':
## 
##     expand

## Loading 'metafor' package (version 2.0-0). For an overview 
## and introduction to the package please type: help(metafor).

source("meta.eval.r") ## Multiple aggregate
  
  
## Simplify the estimates column
meta2 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta2 <- merge(meta2, est, by="Estimate")

 
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.data.csv")
meta2 <- merge(meta2, ms.data, by="uniqueID")
## subset for abundance only
meta2 <- meta2 %>% filter(est.reduced == "diversity")
  
## Select domestic grazers only
meta2 <- meta2 %>% filter(grazer.simplified == "domestic")
 

  
## Create Unique identifier column
meta2[,"UniqueSite"] <- paste(meta2$uniqueID, meta2$Higher.Taxa, meta2$Taxa.1, meta2$Functional.Group, meta2$estimate.simplified, sep="-")

## convert se to sd
meta2[meta2$Stat=="se","Value"] <- meta2[meta2$Stat=="se","Value"] * 1.96
meta2[meta2$Stat=="se","Stat"] <- "sd"

## drop comparisons that are not pairwise
meta2 <-  meta2 %>% filter(grazing.compare.simple == "ungrazed" | grazing.compare.simple == "grazed")
meta2 <- meta2 %>% filter(uniqueID != "30") ## drop study 30 because anomalous 



## Drop plants
meta2 <- meta2 %>% filter(Taxa.1 != "Plants") ## examine animals only


## Use function to extract summary statistics for comparisons
## meta.eval  arguments are (meta.data, compare, ids , stats)
grazed.compare <- meta.eval(meta2, grazing.compare.simple, UniqueSite, Stat)

## Combine the lists into same dataframe
## Rename Columns in second dataframe
grazed.stat <- grazed.compare [[2]] ## extracted statistics 
names(grazed.stat) <- c("UniqueSite","grazed_mean","grazed_sd","ungrazed_mean","ungrazed_sd","grazed_n","ungrazed_n") ## rename columns to match
grazed.raw <- grazed.compare[[1]] ## calculated statistics from raw values

## Join two dataframes
meta.stat <- rbind(grazed.raw, grazed.stat[, names(grazed.raw)])


meta.ready <- escalc(n1i = ungrazed_n, n2i = grazed_n, m1i = ungrazed_mean, m2i = grazed_mean, sd1i = ungrazed_sd, sd2i = grazed_sd, data = meta.stat, measure = "SMD", append = TRUE)

## Warning in ifelse(mi <= 1, NA, exp(lgamma(mi/2) - log(sqrt(mi/2)) -
## lgamma((mi - : value out of range in 'lgamma'

## Warning in ifelse(mi <= 1, NA, exp(lgamma(mi/2) - log(sqrt(mi/2)) -
## lgamma((mi - : value out of range in 'lgamma'

## clean up meta.ready
meta.ready <- na.omit(meta.ready) ## drop NAs
meta.ready <- meta.ready[meta.ready$grazed_mean<80,] ## Drop extreme high variance


## separate out the identifiers
siteID <- matrix(unlist(strsplit(meta.ready$UniqueSite,"-")),ncol=5, byrow=TRUE)
siteID <- data.frame(siteID) ## recreate as dataframe
colnames(siteID) <- c("Study","higher.taxa","taxa","FG","measure") ## add column names
meta.ready <- cbind(data.frame(meta.ready), siteID)

#random-effects meta-analysis for grazed vs ungrazed plots
m1 <- rma(yi=yi, vi=vi,  data = meta.ready)
summary(m1) 

## 
## Random-Effects Model (k = 39; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -67.6191  135.2382  139.2382  142.5134  139.5811  
## 
## tau^2 (estimated amount of total heterogeneity): 1.5561 (SE = 0.4338)
## tau (square root of estimated tau^2 value):      1.2474
## I^2 (total heterogeneity / total variability):   87.24%
## H^2 (total variability / sampling variability):  7.84
## 
## Test for Heterogeneity: 
## Q(df = 38) = 219.2507, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub    
##   0.6382  0.2209  2.8884  0.0039  0.2051  1.0712  **
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## get last year grazed as covariate
yr.grazed <- meta2 %>% group_by(uniqueID) %>% summarize(yr=max(yr.grazed))
colnames(yr.grazed)[1] <- "Study"
site.yr <- merge(siteID, yr.grazed, by="Study")
meta.ready[,"yr.grazed"] <- site.yr$yr

#mixed-effects meta-analysis for grazed vs ungrazed plots
m2 <- rma(yi=yi, vi=vi, mods=~FG+higher.taxa-1,  data = meta.ready)

## Warning in rma(yi = yi, vi = vi, mods = ~FG + higher.taxa - 1, data =
## meta.ready): Redundant predictors dropped from the model.

summary(m2) 

## 
## Mixed-Effects Model (k = 39; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -57.6547  115.3094  129.3094  139.7849  133.7894  
## 
## tau^2 (estimated amount of residual heterogeneity):     1.4519 (SE = 0.4440)
## tau (square root of estimated tau^2 value):             1.2049
## I^2 (residual heterogeneity / unaccounted variability): 84.62%
## H^2 (unaccounted variability / sampling variability):   6.50
## 
## Test for Residual Heterogeneity: 
## QE(df = 33) = 186.9464, p-val < .0001
## 
## Test of Moderators (coefficient(s) 1:6): 
## QM(df = 6) = 16.9346, p-val = 0.0095
## 
## Model Results:
## 
##                        estimate      se     zval    pval    ci.lb   ci.ub
## FGAll                    0.7364  0.3279   2.2460  0.0247   0.0938  1.3789
## FGDetrivorous            0.3734  1.2370   0.3019  0.7627  -2.0510  2.7979
## FGHerbivorous            1.0664  0.5973   1.7853  0.0742  -0.1043  2.2370
## FGPollinator             1.6756  0.5336   3.1400  0.0017   0.6297  2.7214
## FGPredator               0.4683  0.9399   0.4982  0.6184  -1.3740  2.3105
## higher.taxaVertebrate   -0.8916  0.4937  -1.8059  0.0709  -1.8594  0.0761
##                          
## FGAll                   *
## FGDetrivorous            
## FGHerbivorous           .
## FGPollinator           **
## FGPredator               
## higher.taxaVertebrate   .
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## compare inclusion of moderators
(.9539-.9352)/.9352 ## explains an additional 2%

## [1] 0.01999572

## Produce a forest plot to determine the effect sizes for each study
forest(m2)

regtest(m2)

## 
## Regression Test for Funnel Plot Asymmetry
## 
## model:     mixed-effects meta-regression model
## predictor: standard error
## 
## test for funnel plot asymmetry: z = 0.2460, p = 0.8057

confint(m1)

## 
##        estimate   ci.lb   ci.ub
## tau^2    1.5561  0.9783  3.1496
## tau      1.2474  0.9891  1.7747
## I^2(%)  87.2371 81.1213 93.2589
## H^2      7.8352  5.2970 14.8344

## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m2)

## Calculate rosenthals Failsafe number
fsn(yi, vi, data=meta.ready)

## 
## Fail-safe N Calculation Using the Rosenthal Approach 
## 
## Observed Significance Level: <.0001 
## Target Significance Level:   0.05 
## 
## Fail-safe N: 847

# ## generate plot with spaces inbetween
# forest(m1, atransf=exp, cex=0.75, ylim=c(-1, 24),
#        order=order(meta.ready$GI.compare,meta.ready$taxa), rows=c(3:4,7,10:16,19:21),
# #         mlab="RE model for all studies", psize=1, slab= paste(meta.ready$Study, meta.ready$taxa, meta.ready$measure))
# 
# addpoly(res.w, row=18, cex=0.75, atransf=exp, mlab="RE model for green wall")
# addpoly(res.r, row= 9, cex=0.75, atransf=exp, mlab="RE model for green roof")
# addpoly(res.rd, row= 6, cex=0.75, atransf=exp, mlab="RE model for roadsides")
# addpoly(res.p, row= 2, cex=0.75, atransf=exp, mlab="RE model for retention ponds")

Compare regional habitats

meta3 <- meta
est <- read.csv("data//Unique.Estimates.Column.csv")
meta3 <- merge(meta3, est, by="Estimate")

 
## join manuscript data
ms.data <- read.csv("data//synthesisdata//manuscript.data.csv")
meta3 <- merge(meta3, ms.data, by="uniqueID")

library(raster)

## Warning: package 'raster' was built under R version 3.5.2

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.5.2

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

## Convert lat lon to spatial points
meta3 <- meta3 %>% filter(lat != "") ## remove blanks for lat and lon
meta3 <- meta3[!grepl(";", meta3$lat),] ## remove multiple entries
meta3$lat <- as.numeric(levels(meta3$lat))[meta3$lat] ## convert from factor to number

## Warning: NAs introduced by coercion

meta3$lon <- as.numeric(levels(meta3$lon))[meta3$lon] ## convert from factor to number

## Warning: NAs introduced by coercion

## generate dataframe
gps <- data.frame(uniqueID = meta3$uniqueID, lon=meta3$lon, lat=meta3$lat)
gps <- gps[!duplicated(gps$uniqueID),]


## assign spatial points
crs.world <-CRS("+proj=longlat +datum=WGS84")
coordinates(gps) <- ~lon+lat
proj4string(gps) <- crs.world

## load rasters
temp <- raster("C:\\Users\\Fitz\\Downloads\\wc2.0_30s_bio\\bio1.tif")
prec <- raster("C:\\Users\\Fitz\\Downloads\\wc2.0_30s_bio\\bio12.tif")

clim <- stack(temp,prec)


## extract temperature and preciptiation
gps.clim <- data.frame(extract(clim, gps))
gps.clim["Study"] <- gps$uniqueID


## compare climate & grazing duration variables against effect size

## get last year grazed as covariate
yr.grazed <- meta2 %>% group_by(uniqueID) %>% summarize(yr=max(yr.grazed))
colnames(yr.grazed)[1] <- "Study"
site.yr <- merge(siteID, yr.grazed, by="Study")
meta.ready[,"yr.grazed"] <- site.yr$yr

#mixed-effects meta-analysis for grazed vs ungrazed plots
m3 <- rma(yi=yi, vi=vi, mods=~yr.grazed,  data = meta.ready)

## Warning in rma(yi = yi, vi = vi, mods = ~yr.grazed, data = meta.ready):
## Studies with NAs omitted from model fitting.

summary(m3) 

## 
## Mixed-Effects Model (k = 28; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -47.6117   95.2234  101.2234  104.9977  102.3143  
## 
## tau^2 (estimated amount of residual heterogeneity):     1.8660 (SE = 0.6103)
## tau (square root of estimated tau^2 value):             1.3660
## I^2 (residual heterogeneity / unaccounted variability): 87.58%
## H^2 (unaccounted variability / sampling variability):   8.05
## R^2 (amount of heterogeneity accounted for):            4.41%
## 
## Test for Residual Heterogeneity: 
## QE(df = 26) = 180.9830, p-val < .0001
## 
## Test of Moderators (coefficient(s) 2): 
## QM(df = 1) = 2.1745, p-val = 0.1403
## 
## Model Results:
## 
##            estimate      se     zval    pval    ci.lb   ci.ub    
## intrcpt      1.1877  0.4102   2.8956  0.0038   0.3838  1.9916  **
## yr.grazed   -0.0262  0.0178  -1.4746  0.1403  -0.0611  0.0086    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## merge climate with dataset
site.clim <- merge(siteID, gps.clim, by.x="Study", all.x=TRUE)
meta.ready[,"temp"] <- site.clim$bio1
meta.ready[,"prec"] <- site.clim$bio12

#mixed-effects meta-analysis for grazed vs ungrazed plots
m4 <- rma(yi=yi, vi=vi, mods=~temp+prec-1,  data = meta.ready)

## Warning in rma(yi = yi, vi = vi, mods = ~temp + prec - 1, data =
## meta.ready): Studies with NAs omitted from model fitting.

summary(m4) 

## 
## Mixed-Effects Model (k = 37; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -64.3611  128.7222  134.7222  139.3883  135.4964  
## 
## tau^2 (estimated amount of residual heterogeneity):     1.8067 (SE = 0.5139)
## tau (square root of estimated tau^2 value):             1.3441
## I^2 (residual heterogeneity / unaccounted variability): 88.80%
## H^2 (unaccounted variability / sampling variability):   8.93
## 
## Test for Residual Heterogeneity: 
## QE(df = 35) = 232.2490, p-val < .0001
## 
## Test of Moderators (coefficient(s) 1:2): 
## QM(df = 2) = 6.2829, p-val = 0.0432
## 
## Model Results:
## 
##       estimate      se    zval    pval    ci.lb   ci.ub   
## temp    0.0225  0.0391  0.5756  0.5649  -0.0542  0.0992   
## prec    0.0006  0.0005  1.0012  0.3168  -0.0005  0.0016   
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1